The sign problem in Monte Carlo simulations of frustrated quantum spin systems 
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We discuss the sign problem arising in Monte Carlo simulations of frustrated quantum spin sys- 
tems. We show that for a class of "semi-frustrated" systems (Heisenberg models with ferromagnetic 
couplings Jz(r) < along the z-axis and antiferromagnetic couplings Jxy(r) — —Jz{r) in the xy- 
plane, for arbitrary distances r) the sign problem present for algorithms operating in the z-basis 
can be solved within a recent "operator-loop" formulation of the stochastic series expansion method 
(a cluster algorithm for sampling the diagonal matrix elements of the power series expansion of 
exp(— /3iif) to all orders). The solution relies on identification of operator-loops which change the 
configuration sign when updated ( "merons" ) and is similar to the meron-cluster algorithm recently 
proposed by Chandrasekharan and Wiese for solving the sign problem for a class of fermion models 
(Phys. Rev. Lett. 83, 3116 (1999)). Some important expectation values, e.g., the internal energy, 
can be evaluated in the subspace with no merons, where the weight function is positive definite. 
Calculations of other expectation values require sampling of configurations with only a small number 
of merons (typically zero or two), with an accompanying sign problem which is not serious. We also 
discuss problems which arise in applying the meron concept to more general quantum spin models 
with frustrated interactions. 

PACS numbers; 75.40.Mg, 75.10.Jm, 02.70.Lq 
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I. INTRODUCTION 

Recently, there have been several significant develop- 
ments of more efficient Monte Carlamethods for interact- 
ing quantum maiw-body systems Jj The Trotter decom- 
position formulafl'B has traditionally been used as a start- 
ing point for finite-tepiperature simulation algorithms, 
such as the worldlinda and fermion determinantu meth- 
ods. It introduces a systematic error that can be removed 
only by carrying out simulations for several different 
imaginary time discretizations At and subsequently ex- 
trapolating to At — 0. Such extrapolations are not nec^ 
essary with the stochastic series expansion method,ErLl 
which is based on sampling the power series expansion of 
exp{—(3H) to all orders and is related to a less^ gejieral 
method proposed much earlier by Handscomb.l3il3 Re- 
sults that are exact to within statistical-erEors can also be 
directly obtaiaed with recent worldlinetll^Ej and fermion 
determinantlij algorithms formulated in continuous imag- 
inary time. Even more significant than the elimination 
of the Trotter decnmppsition error are generalizations to 
the quantum caseli3ll3'Ej'LLof cluster algorithms developed 
for classical Monte Carlo.liZl These "loop algorithms" (so 
called because the clusters are loops on a space-time lat- 
tice) can reduce the autocorrelation times by orders of 
magnitude and enable highly accurate studies of systems 
in parameter regimes where previous algorithms encoun- 
tered difficulties due to long autocorrelation and equili- 
bration times. 

In spite of these developments, the class of mod- 
els which can be studied using quantum Monte Carlo 
methods is, -still severely restricted due to the "sign 
problem" Jlao i.e., the non-positive-definiteness of the 
weight function that can arise in transforming a quan- 



tum problem into a form resembling a classical statisti- 
cal mechanics problem. There are two classes of systems 
for which this issue is particularly pressing — interacting 
fermions in more than one dimension and quantum spin 
systems with frustrated interactions (in any number of 
dimensions). For fermions in one dimension, and hop- 
ping between nearest-neighbor sites only, the sign prob- 
lem can be avoided because the fermion anticommutation 
relations do not come into play (other than introducing 
a hard-core constraint) in the one-dimensional real-space 
path integral. In two or more dimensions (or even in 
one dimension if hopping further than between nearest 
neighbors is included), permutations of fermions during 
the propagation in imaginary time lead to a mixed-sign 
path integral which typically cannot be efficiently evalu- 
ated using Monte Carlo methods. The sign problem can 
be avoided with the fermion determinant algorithm in 
special cases, such as the half-filled Hubbard model (be- 
cause of particle-hole symmetry) jj but in other cases sim- 
ulations are i^stricted to high temperatures and/or small 



system sizes£3 For frustrated spin systems the source of 
the sign problem is different. A minus sign appears for 
every event in the path integral in whichj-two antiferro- 
magnetically interacting spins are flipped.E3 This causes 
an over-all minus sign if the total number of spin flips is 
odd, which can be the case, e.g., for a triangular lattice 
or a square lattice with both nearest- and next-nearest- 
neighbor interactions. Simulations of quantum spin sys- 
tems are therefore restricted to models with no frustra- 
tion (in the off-diagonal part of the Hamiltonian), such as 
ferromagnets, or antiferromagnets on bipartite lattices. 

A promising approach to solving the sign problem 
was recently suggested by Chandrasekharan and Wiese.E2l 
They considered a system of spinless fermions on a two- 
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dimensional square lattice, within the context of the 
worldline loop algorithm.Ej They showed that, for this 
particular model and for a certain range of nearest- 
neighbor repulsion strengths, the properties of the loops 
can be used to eliminate the sign problem. Flipping a 
loop can change the number of fermion permutations 
from odd to even, or vice versa, thereby also changing 
the over-all sign of the configuration. Such sign-changing 
loops are called "merons" . The magnitude of the weight 
is not affected by flipping a meron and therefore all con- 
figurations with one or more merons cancel in the par- 
tition function. The subspace of zero merons is posi- 
tive definite and can be sampled without a sign problem. 
Typical operator expectation values of interest also con- 
tain contributions from configurations with two merons 
which therefore also have to be included in the simu- 
lation and introduce a "mild" sign problem. The rela- 
tive weights of the zero- and two-meron subspaces to be 
sampled can further be chosen in an optimum way using 
a reweighting technique, which further reduces the sign 
problem. 

In this paper we explore an analogous method for solv- 
ing the sign problem for frustrated quantum spin models. 
We consider the operator-looe formulation of the stochas- 
tic series expansion method,Ll in which sequences of two- 
spin operators are sampled by forming clusters (loops) of 
operators that can be simultaneously updated without 
changes in the weight function. The updated clusters 
contain operators acting on the same spins, but diago- 
nal operators may be changed to off-diagonal ones and 
vice versa. For a model with frustrated interactions an 
operator-loop update can lead to a sigp-change. In anal- 
ogy with Chandrasekharan and Wiesec3 we will refer to 
such sign-changing loops as "merons" . The sign problem 
can be solved if the operator-loops for a given configura- 
tion can be uniquely defined and the weight function is 
positive definite in the configuration subspace containing 
no merons. Unfortunately, we find that these criteria are 
in general difficult to satisfy. Operator-loop algorithms 
with uniquely determined loops are typically non-ergodic 
for frustrated systems, and with supplemental local up- 
dates for ergodicity there are mixed signs in the zero- 
meron subspace. In fact, in such cases merons typically 
do not even exist, i.e., none of the operator-loops can 
change the sign when flipped. We have found only one 
spin system for which the sign problem can be eliminated 
using merons — the Heisenberg model with ferromag- 
netic couplings Jz (r) < along the z-axis and frustrated 
antiferromagnetic couplings Jxyif) = —Jz{t) in the plane 
perpendicular to this axis, i.e., the Hamiltonian 

^ = - E - l^^tSj + (1) 

where > and the range of the couplings is arbi- 
trary. We have implemented a meron algorithm for this 
model on a square lattice with nearest- and next- nearest- 
neighbor couplings J(l) and J(\/2). Standard algorithms 



for this model have a severe sign problem when using the 
z-axis as the quantization axis, however, it can be avoided 
by a simple rotation to the ^-representation. Using the 
SSE algorithm and the meron concept, the sign prob- 
lem can be eliminated also in the z representation. With 
both representations accessible in simulations, correla- 
tion functions both parallel and perpendicular to the z 
direction can be easily evaluated. 

The model, Eq. (|l|), can be mapped onto a hard-core 
boson model with attractive interactions and frustrated 
hopping. Frustration in the potential energy has been in- 
vestigated in this context as a possible raachanism to ren- 
der a disordered bosonic ground statcEil Frustration in 
the hopping [the xy-term in Eq. (|l|)] should decrease the 
tendency to forming off-diagonal long-range order and 
could then lead to a normal fluid (non-superfluid) . How- 
ever, the highly symmetric case considered here has a 
trivial, ordered ground state; the fully polarized ferro- 
magnetic state (corresponding to a completely filled lat- 
tice of hard-core bosons; a trivial case of normal solid). 
Effects of frustration only come into play at finite temper- 
ature, where the model is different from the correspond- 
ing isotropic ferromagnet (on non frustrated, bipartite 
lattices the two models are equivalent, since the sign of 
the xy-teim. can be switched by a spin rotation on one of 
the sublattices) . 

Although we have not been able to solve the sign 
problem for other cases, such as the Heisenberg model 
with completely antiferromagnetic interactions [Jz(l) = 
Jxy{^) > and Jz{V^) = JxyiV^) > 0], our work never- 
theless gives some further insights into the meron concept 
and what is required in order to solve the sign problem 
for arbitrary couplings. 

The outline of the rest of the paper is the following: In 
Sec. II we review the basics of the stochastic series expan- 
sion method and discuss operator-loop updating schemes 
for both ferromagnetic and antiferromagnetic couplings. 
In Sec. Ill we present the solution of the sign problem for 
the Jz{r) = —Jxy{r) model. The reweighting technique 
is analyzed in some detail in Sec. IV. In Sec. V we discuss 
some simulation results for the semifrustrated model and 
make comparisons with the isotropic Heisenberg ferro- 
magnet. We summarize our work in Sec. VI. 

II. OPERATOR-LOOP ALGORITHM 

In this section we first briefly review the expansion un- 
derlying the SSE method and then discuss the operator- 
loop updates used to efficiently sample the expansion. 
We here assume a non-frustrated case and postpone the 
discussion of the sign problem for frustrated models to 
Sec. HI. For definiteness we consider the S = 1/2 Heisen- 
berg model 

H = ±jY,S^-S,, (2) 
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where (ij) denotes a pair of nearest-neighbor spins on 
a cubic lattice (in an arbitrary number of dimensions), 
and J > 0. Depending on the sign, the model is an 
antifcrromagnet (— ) or a ferromagnet (+). To construct 
the SSE configuration space the Hamiltonian is rewritten 
as a sum of diagonal and off-diagonal operators 



J 



M 



(3) 
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where the index b denotes an interacting spin pair (bond) 
{i{b),j{b)) and C is an irrelevant constant equal to MJ/A, 
where M is the total number of pairs of interacting spins. 
The bond-indexed operators are given by 



(4) 
(5) 



Note that the eigenvalues of both the diagonal (i?i,b) and 
the off-diagonal (i?2,b) operators are and 1, both for 
the antifcrromagnet and the ferromagnet. Tie partition 
function Z = Tr{exp(— /3i/)} is expanded as 



(6) 



Q n— 



in the basis {\a)} = {\Sl,S^, . . . ,S^)}, where N is the 
number of spins. Terms of order greater than n ^ N [3 
give an exponentially vanishing contribution and for the 
purpose of a stochastic sampling the expansion can there- 
fore be truncated at some n — L of this order without 
loss of accuracy (see, e.g., Rcf. ^ for details on how to 
choose a sufficiently high truncation power). Additional 
unit operators Hq.o = 1 are introduced to rewrite Eq. (|^) 
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where Sl denotes a sequence of operator-indices: 



Sl = (ai,6i)i, (02,62)2 



iaL,b 



L)L, 



(7) 



(8) 



with a^ G {1,2} and b, € {1,...,M}, or {a^,b^) = (0,0). 
The number of non-(0,0) elements in 5*^ is denoted n, 
while n2 denotes the number of off-diagonal operators 
in the sequence. Note that since the expectation value 
in Eq. (|^) is always equal to zero or one, the sign of a 
term is negative only if 712 is odd. This sign problem 
occurs (only) when frustration is present and is the main 
topic of this paper. However, for the discussion of the 
sampling procedures, in this section we assume a positive 
definite expansion. We introduce the notation \a{p)) for 
a propagated state 



Hp)) ^tlHa^M, 



i=l 



where for an allowed configuration |a(0)) ~ 
and the weight function corresponding to (^ 



(9) 



a{L)) = \a) 
) is given by 



(JI3y^(L-n)\ 



(10) 



Having established the framework we will proceed 
to describe the procedures for importance sampling 
of the terms (a, S'l) according to the weight (px|). 
The initial state can be a sequence of the form 
(0,0)1, (0,0)2, ■ ■ • , (0,0)l (subscripts on the index pairs 
will sometimes be used to denote the position in 5*^) 
and a random state \a). An ergodic procedure for sam- 
pling the terms is achieved using two types of basic up- 
dates; a simple substitution of single diagonal operators 
and the operator-loop update which involves simultane- 
ous updates of a number (in principle varying between 1 
and n) of diagonal and off-diagonal operators. 

The diagonal update is carried out by traversing the 
operator- index sequence Sl from beginning {p — 1) 
to end (jp = L). Operator substitutions of the form 
(O, 0)p ^ (1, b)p are attempted where possible, while the 
stored state \a) is updated every time an off-diagonal op- 
erator is encountered so that the state |a(p)) is available 
when needed. With the weight function (10) detailed bal- 
ance can be seen to be satisfied if the acceptance proba- 
bilities are taken to be 

P [(0, 0), ^ (1, 6) J = ^Pi-^iP)\HMp)) ^ (11) 



P[(l,6)p^(0,0)J =- 



L — n 
L-n + 1 



Mp{abip)\Hi^b\abip))' 



(12) 



Note that the diagonal update changes the expansion 
power n by ±1. Off-diagonal operators cannot be in- 
troduced one-by-one because of the periodicity condition 
\a{L)) — |a(0)). Local updaips involving two operators 
can be used for this purpose,El but are more complicated 
and far less efficient than the operator-loop update,Ll 
which is discussed next. 

We use a largely pictorial description of the operator 
loops. First we consider the antifcrromagnet. Note that 
in this case the only non-zero matrix elements of the bond 
operators are 



at |i^iUT) = (Til^ilTi) = i, 

at m u) = ill \H2\ it) = 1, 



(13) 



i.e., they can act only on antiparallel spins. An example 
of a term in the expansion for a four-site antifcrromagnet 
is depicted in Fig. (0) . This representation makes evident 
the close relationship between the SSE expansion and the 
Euclidean path integral. An imaginary time separation t 
corresponds to a distribution of propagatijMis Ap between 
states, centered around Ap = {t / 0)nn't3 We will for 
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|a(0)> A_ 
|a(l)> A 
|a(2)> A 



A 
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|a(3)> A V V 
|a(4)> A V 
|a(0)> A V 



A 
A 



A V 



Hi.i 

H2,3 

Hi.i 

H2,3 
Hl,2 



FIG. 1. Representation of a term in the SSE expansion of 
a four-site antiferromagnet. Up and down spins are repre- 
sented as solid and open triangles, respectively. The horizon- 
tal bars indicate the presence of diagonal (narrow lines) and 
off-diagonal (thick lines) operators. 



convenience here refer to the propagation as the time 
dimension. |— ■ 

The general ideate behind the loop update is to flip a 
cluster of spins in the configuration in such a way that the 
weight is not changed. With the SSE method there will 
also have to be changes made to the operators acting on 
the spins, since otherwise operators Hi or H2 may act on 
parallel spins, resulting in zero- valued matrix elements. 
Since one of the states \a{p)) and the operator sequence 
Sl uniquely define the whole spin configuration, the SSE 
loops are in practice treated as loops of operators, the 
exact meaning of which will be made clear below. 



▲ V 



A V 



FIG. 2. Spin flip and accompanying operator exchange dur- 
ing a loop update for an antiferromagnet. The dashed line 
indicates part of the loop. 

Consider one of the operators Hi i in Fig. (|l|). It can be 
depicted as a "vertex" with four legs associated with spin 
states I or "f . If we flip the upper left spin, a vanishing 
matrix element results. But if we flip both upper (or 
lower) spins and also change the operator type to off- 
diagonal H2.1, an allowed matrix element is generated; 
see Fig. (||) . Using this idea we can form a cluster of spins 
by choosing a random spin S^ (p) in the configuration and 
traversing up or down until we encounter an operator 
(bond) acting on that spin, then switch to the second 
spin of the bond and change the direction of traversing 
the list. Eventually we will necessarily arrive back at 
the initial starting point, whereupon a closed loop has 
formed. All the spins on this loop can be flipped if the 
operators encountered are also switched [(1,6) <-> (2,6)]. 
Note that the same operator can be encountered twice, 
which results in no net change of operator type (but the 
spins at all four vertex legs are flipped). 

The whole configuration can be uniquely divided up 
into a set of loops, so that each spin belongs to one and 
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FIG. 3. The SSE space-time configuration of Fig. 1 
uniquely divided up into loops (left). The right hand configu- 
ration results from flipping the loop indicated by the dashed 
line. Note the periodic boundary conditions in the vertical 
("time") direction. 



only one cluster; see Fig. (|3|), where our example con- 
figuration has been divided up into three clusters. All 
loops can be fiipped independently with probability 1/2 
— in Fig. (||) we depict the result of hipping the largest 
loop of the example configuration. A full operator-loop 
update amounts to dividing up the configuration into 
all of its loops which arc flipped with probability 1/2. 
The (random) decision of whether or not to flip a loop 
can be made before the loop is constructed, so that each 
loop has to be traversed only once. Spin "lines" Sf{p), 
p = 0, . . . , n, which are not acted upon by any of the op- 
erators in Sl will not be included in any of the operator- 
loops. They correspond to "free" spins which can be 
flipped with probability 1/2. Such a line can also be con- 
sidered a loop, and then it will always be true that every 
spin S^{p) belongs to one loop. Free spins appear fre- 
quently at high temperatures, when the total number of 
operators n ^ TV, but are rare at low temperatures. 

Note that the spin states at the four legs of the 
operator-vertices completely determine the full spin con- 
figuration, except for free spins that happen not to be 
acted upon by any of the operators in Sl- Hence, the 
operator-loop update can be carried out using only a 
linked list of of operators, i.e., an array of vertices with 
four spin states and associated pointers to the "previ- 
ous" and "next" vertex associated with the same spin. 
The storage requirements and the number of operations 
needed for carrying out a full operator-loop update then 
scale as ~ N(3 instead of ~ N'^f3 if the full spin configu- 
ration were to be used. 

In a simulation we first make a full cycle of diagonal 
updates in the sequence Sl and then create the linked 
list of vertices in which the operator-loop updates are 
carried out. The vertex list is then mapped back onto 
the sequence Sl and the state |a(0)). Alternatively, the 
linked list can be updated simultaneously with each di- 
agonal operator substitution, so that it does not have to 
be re-created for each Monte Carlo step — depending on 
the model studied there may be significant differences in 
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execution time between the two approaches. 

▲ ▲ V ▲ 



▲ A 



▲ V 



FIG. 4. Spin flip and accompanying operator exchange dur- 
ing a loop update for a ferromagnet. The dashed line indicates 
part of the loop. 




FIG. 5. SSE space-time configurations uniquely divided 
into loops for the case of a ferromagnet. The left and right 
hand configurations differ by fiipping the loop indicated by a 
dashed line. 

For the ferromagnet we can construct the loops in a 
similar manner, but the non-zero matrix elements are 
now 



(ii li^il 

(it 1^2 I 



ii) 

Ti> 



(TT \Hi\ 

(Ti 1^2 I 



TT> 
it) 



(14) 



i.e., the off-diagonal operators act only on antiparal- 
lel spins, as before, whereas the diagonal ones can act 
only on parallel configurations. This implies qualita- 
tive changes in the structure of the loops, as depicted 
in Fig. (^. If we again consider an operator Hi and flip 
the upper left spin, we note that we need to flip the lower 
right spin and change the operator to H2. Hence, instead 
of changing the direction of traversing the configuration 
every time an operator is encountered we now continue 
in the same direction after switching to the second spin 
of the bond. Any configuration can still be uniquely di- 
vided up into loops that can be flipped with probability 
1/2. An example of a ferromagnetic configuration with 
its loops is shown in Fig. (||). 

Note that since the loops for the ferromagnet never 
change direction as they go through the lattice, every 
single loop has to traverse each state \a{p)) at least once. 
It follows that the number of sites N is an upper bound 
of the number of loops. The antiferromagnetic loops, on 
the other hand, traverse the lattice in both directions and 
the number of loops is therefore instead limited by the 
total number of operators n ^ Nf3. As a consequence 
of the change of direction, for the antiferromagnet the 



linked list of vertices must to be bi-directional, but for 
the ferromagnet it is sufficient to keep a singly-directional 
list. 

The diagonal and operator-loop updates satisfy de- 
tailed balance and the combination of them leads to er- 
godic sampling for a ferromagnet on any lattice, and for 
antiferromagnets on bipartite lattices — for frustrated 
antiferromagnets there are complications, in addition to 
the sign problem, which will be discussed further in the 
next section. The operator-loop sampling is highly ef- 
ficient, with integrated auto-correlation times typically 
less than one Monte Carlo step. 

The loop construction described here relies on the rota- 
tional invariance of the models, i.e., the fact that both the 
diagonal and off-diagonal matrix elements in Eqs. (|l^ ) 
and (|l^) are equal to one. For a general anisotropic case 
the loops will lead to weight changes when flipped and 
must then be assigned "a-posteriori" acceptance prob- 
abilities which typically become small for large lattices 
at low temperatures. Other Ispes of loops avoiding this 
problem have been proposed,!] but will not be discussed 
here. 



III. THE SIGN PROBLEM 

The notorious sign problem arises in stochastic sam- 
pling when the function used to weight the different con- 
figurations in not positive definite. A typical quantity 
that can be calculated by Monte Carlo (importance sam- 
pling) is of the form 



Y.^A{x,)W{x,) 



{A{x)) 



(15) 



where W is the weight function and A the measured 
quantity, which both depend on the general coordinate 
X of the configuration space sampled. When the coordi- 
nates are sampled according to relative weight, the de- 
sired quantity is simply the arithmetic average of the 
measured function A{x)^ as indicated by the notation 
{A{x))w above. If the weight function is not positive 
definite, the sampling can be done using the absolute 
value of the weight, and the expectation value can be 
calculated according to 



{A{x)s{x))m\ 



,{x)) 



w\ 



(16) 



where s{x) equals ±1, depending on whether the sign 
of the weight function is positive or negative. For most 
models, where a sign problem is present, the average sign 
{s(x))\\Y\ decreases exponentially to zero as a function 
of inverse temperature and system size, and the relative 
statistical errors of calculated quantities increase expo- 
nentially. 

A meron-cluster solution to the sign problem using 
loop updates of fermion world-line configurations ^¥as re- 
cently proposed by Chandrasekharan and Wiese.Ey This 
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approach is based on the idea that if it is possible to 
map every configuration with negative weight uniquely 
to a corresponding configuration with equal weight mag- 
nitude but opposite sign, then the partition function can 
be sampled without a sign problem simply by not in- 
cluding any configuration which is a member of such a 
canceling pair. In the meron-cluster algorithm, flipping a 
loop of spins can lead to a sign change without change in 
the magnitude of the weight, and such a "meron" hence 
identifies a canceling pair of configurations. Here we will 
present a similar approach within the SSE operator-loop 
method for frustrated quantum spins. 

Let us consider the Heisenberg antiferromagnet dis- 
cussed in the previous section. From Eq. (Q) we see that 
a configuration has negative weight if the total number 712 
of off-diagonal operators is odd. This can only be the case 
on frustrated lattices. As described in the previous sec- 
tion, any SSE configuration can be divided up uniquely 
into a number of loops. Flipping a loop interchanges the 
diagonal and off-diagonal operators, but leaves the total 
number of operators unchanged. It follows that the sign 
will change if and only if a loop passing through an odd 
number of operators is flipped (two passes through the 
same operator is counted as two operators). Since the 
total weight remains unchanged we have thus found the 
desired mapping between positive and negative config- 
urations (assuming that there exist loops which change 
the sign when flipped, which in fact is not always the 
case) . In analogy with previous work we call such a sign- 
flipping loop a "meron" . Let us now see how we can use 
this concept to measure observables. 

As in Ref. we consider improved estimators that 
average over all loop configurations. Denote the number 
of loops in the system Nl- Since each loop can be in one 
of two states there is a total of 2^^ configurations that 
can be reached by flipping all combinations of the loops 
present. The improved estimate therefore takes the form 



{A) = 



{{A{x)s{x)))\w\ 
{{s{x)))\w\ 



(17) 



where the double brackets denote an average over all the 
loop states for each generated SSE configuration, e.g., 



(18) 



The general coordinate x here refers to the SSE configu- 
ration space (a, Sl) and xi refers to one out of the 2^^ 
possible outcomes of "flipping" a number of loops. 

Let us consider this average. Denote the state of a 
loop with 5, with two possible "orientations" 5 £ {"f, J,}. 
Since flipping one loop does not affect any other loops (in 
terms of their paths taken), the sign of a configuration 
factors according to 



Nl 



3(61,62, ...,6nl) = Y[s{6,), 



(19) 



where s{6) = ±s{5), where 6 denotes a flipped loop, and 
the sign is negative for merons and positive otherwise. 
Since flipping any meron leads to two terms that cancel, 
it follows that 



2^ *(^') = 
1=1 



(20) 



where hm denotes the total number of morons. The sign 
in front of the delta function is the "inherent" sign of the 
configuration, independent of the loop orientation when 
there are no merons present. This sign has to be positive 
for the meron solution to be applicable in practice, and 
then the partition function can be sampled in the positive 
definite subspace of configurations with no merons. 

Having found an expression for the denominator in 
Eq. we need to consider the numerator for cases of 
interest. SSE estimators for a number of important op- 
erators have been discussed, e.g., in Ref. ^. The internal 
energy is given by 



E 



'n)w, 



(21) 



where n denotes the total number of operators in the 
sequence Sl. This number is not affected by the loop 
updates, and hence it follows that 

1 {{n{x)s{x)))\w\ ^ 1 {n{x)6nM,o)\w\ ^^2) 



f3 {{s{x))),w 



P {6nM,o}\W\ 



Assuming that this sector has positive definite weight we 
have therefore completely eliminated the sign problem by 
restricting the simulation to the zero-meron sector. The 
energy is then simply given by 



E = -- 



lix')) 



w 



(23) 



where the superscript indicates the restriction of the 
simulation to the zero-meron sector. 

Next we will consider the magnetic susceptibility. 




/N = I3{M^)/N. 



(24) 



Since M is conserved by the Hamiltonian its value is 
the same in all propagated states; M{p) — Af(0) = M 
{p = I, . . . ,n). In a configuration uniquely divided up 
into loops, every spin Sf{p) belongs to one and only one 
loop, if we count as loops also all "lines" of free spins, 
i.e., the spins Sf(p),p=l,...,n for all sites i which are 
not associated with any operator in the sequence (and 
therefore can be flipped). It follows that the change in 
M{p) when flipping a loop must be the same for all p, 
and hence only loops that go through all states \a{p)) (at 
least once) can change M when flipped. We can there- 
fore introduce a loop magnetization ttil, which is simply 
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equal to the sum of the spins traversed by the loop for 
an arbitrary \a{p)). In the estimator ([l^ ) corresponding 
to the susceptibility the numerator can hence be written 
as 



\ 1=1 



(25) 



The magnetization on a loop always changes sign when 
a loop is flipped; the over-all sign s{xi) only changes sign 
when a meron is flipped. Therefore, in summing over all 
loops in (|2^), a non-zero value results only if the conflg- 
uration has zero or two merons. The full susceptibility 
estimator therefore takes the form 



X 



nM;2 



M,0/ 



(26) 



where Mi and M2 are the indices of the loops correspond- 
ing to merons in a two-meron conflguration. Hence, un- 
like in the case of the total energy, the sign problem has 
here not been completely eliminated, since the zero- and 
two-meron configuration contribute 1 and 0, respectively, 
to the average sign. When the SSE conflguration volume 
grows the relative weight of the zero-meron sector should 
diminish, leading to a decreasing average sign. Chan- 
drasekharan and Wiese stated that the statistical fluctu- 
ation in the improved susceptibility estimator increases 
quadratically with NP^ i.e., much slower than the conven- 
tional exponential increase. They also pointed out that 
this remaining sign problem can be solved by reweight- 
ing the zero- and two-meron sectors with external weight 
factors 'w{0) and w(2). This changes the above formula 
to 

J2i=i \'^i\^nM,ow{2) -I- 2|mAfi||mM2l<5„„,2U'(0) 



('5»M.ou;(2)) 



(27) 



In the next section we will say more about reweighting. 

As we have shown above, the meron concept within 
the SSE method formally leads to exactly the same equa- 
tions as in the world-line simulations of fermion systems 
considered in Ref. |2^. The difference is only in the struc- 
ture of the meron itself; the fermionic meron changes the 
number of particle permutations from even to odd, or 
vice versa, whereas the SSE meron in the case of a frus- 
trated spin system instead changes the number of anti- 
ferromagnetic spin flips from even to odd or vice versa. 
Applying the SSE operator-loop algorithm to a fermion 
system would lead to merons of exactly the same kind 
as those existing within the world-line framework, and, 
conversely, applying a world-line loop algorithm to a frus- 
trated spin system should lead to merons similar to those 
discussed here (there are no diagonal operators in the 
world-line conflgurations, but spin flip events correspond 
to the SSE off-diagonal operators and can change from 



even to odd, or vice versa, in loop updates) . These simi- 
larities are not surprising, considering the close relation- 
ship betateen the SSE expansion and the Euclidean path 
integral .Ej 



FIG. 6. Three off-diagonal operators acting on a triangle 
of spins. The left and right hand configurations differ by 
flipping the loop indicated by a dashed line. Note the periodic 
boundary conditions. 

Now consider the application of the above ideas to the 
Heisenberg model on a square lattice with nearest- and 
next-nearest-neighbor coupHngs J(l) > and J(a/2) > 
(antiferromagnetic). This model has a sign problem since 
the total number 712 of spin-flipping operators in a con- 
flguration can be odd, e.g., three operators on a triangle 
of spins, as shown in Fig. (^. Already this simple exam- 
ple illustrates that the operator-loop algorithm discussed 
above does not sample the full conflguration space and 
that the meron concept therefore cannot be applied to 
solving the sign problem. Since all three bonds are anti- 
ferromagnetic, a loop will change direction every time an 
operator is encountered. In order for the loop to close, 
it therefore has to pass through an even number of op- 
erators and hence flipping the loop cannot change the 
number rt2 of off-diagonal operators from odd to even, or 
vice versa. This is illustrated in Fig. (^, where the only 
effect of flipping the single loop in the system is to ffip 
all the spins; the operators remain unchanged. Hence, 
merons do not even exist within the operator-loop algo- 
rithm for this model, and the sampling is non-ergodic. A 
local update can in principle be used in combination with 
the operator-loops in order to make the sampling ergodic. 
However, a configuration can then have a negative sign 
(of which Fig. (^ is an example) even though there are no 
merons present. The principal requirement of positive- 
definiteness in the zero-meron subspace (which in this 
case is the full space) is hence not fulfilled. A similar 
problem seems to affect all models with frustration in all 
of the spin components. It is possible that some other 
way of constructing the loops could remedy this, e.g., by 
switching to some other, non-trivial basis in which the 
SSE expansion could be carried out. Other ways pro- 
posed for constructing loops in the standard basis, con- 
sidered here do not uniquely define the set of loopaJ and 
can therefore not easily be used with the meron concept. 

We have found one class of spin models for which the 
meron ideas can be successfully applied to solve the sign 
problem: Heisenberg models which are antiferromagnetic 
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in the xy-plane but ferromagnetic along the z-axis., i.e., 
the "semi-frustrated" model (|^) , which can be written as 



M 



where the bond-indexed operators are given by 



Hi.h 

H2.b 



4 + '^i(6)'^i(6) 



(28) 

(29) 
(30) 



On a non-frustrated lattice this model is equivalent to an 
isotropic Heisenberg ferromagnet, since n2 is always even 
and the sign in front of the operators i?2,& in Eq- (|2^ ) is 
irrelevant as (—1)"^ = 1 in Eq. (|^) — the sign can also 
be transformed away by a spin rotation on one of the 
two sublattices. On a frustrated lattice, on the other 
hand, n2 can be odd (the lattice is no longer bipartite 
so that the transformation mentioned above does not re- 
move all signs), and the system is no longer equivalent 
to the isotropic ferromagnet. The model has a classi- 
cal two-fold degenerate ferromagnetic ground state, but 
at finite temperatures the transverse spin components 
are frustrated, and the behavior will be different from 
the isotropic ferromagnet. When simulated in the z- 
basis using standard algorithms the semifrustrated model 
has a severe sign problem, but the zero-meron sector is 
positive-definite and the SSE meron solution can be ap- 
plied. The structure of the operator-loops is the same as 
for the ferromagnet and the loop algorithm is therefore 
ergodic. 

The meron solution can be implemented in several dif- 
ferent ways and we briefly describe how it was done in 
this work: During the sequential diagonal updates the 
linked list, representing the loop structure, is updated si- 
multaneously with each accepted diagonal update. The 
loops are numbered and information is stored on whether 
each loop is a meron or not. During an attempted diago- 
nal update only the loops directly affected by the opera- 
tor substitution are updated. This permits easy and fast 
checking of whether the number of merons in the system 
has changed or not. If the new number of merons is dif- 
ferent form zero or two the update is rejected, whereas if 
the number of merons changes from i to j it is accepted 
with probabihty w{j)/w{i), where i,j G {0,2}, and w{i) 
is the re- weighting factor assigned to meron sector i. 



A. V 



V 
A A 



FIG. 7. Spin flip and accompanying operator exchange 
during a loop update for model with a ferromagnetic 
a;-component, but antiferromagnetic y and z-components. 
The dashed line indicates part of the loop. 

Note that the sign problem for the semifrustrated 
model can also be very easily transformed away by ro- 



tating the ferromagnetic component to the y-direction. 
The Hamiltonian then takes the form 



M 



b=l 



where the bond-indexed operators are given by 



QZ QZ 



^ib-^yj '-'i(b)'^j{b) 



H2,b - ^m^tib) + ^m^jib}' 



(31) 



(32) 
(33) 



and the fundamental spin-flips and operator exchange 
during a loop update is shown in Fig. (0). Being able to 
work in both bases we can easily measure all components 
of the susceptibility. Our main motivation for studying 
this model is to illustrate how the sign problem can be 
removed in the z-basis. Nevertheless, we will also show 
some results calculated in the x-basis. 



IV. REWEIGHTING 

An important technical aspect of the meron-solution is 
the reweighting of the zero- and two-meron sectors, which 
was briefly mentioned in the previous section. Eq. ( [27| ) 
gives the correct estimator for the susceptibility after 
reweighting, but it gives no information on how to do 
the reweighting in practice. How to determine the op- 
timal reweighting and whether reweighting changes the 
scaling of the relative error are important questions to be 
considered in this section. 



- 

50 100 150 200 

MC steps 

FIG. 8. Fluctuations of n during a random process with 
two outcomes (n — 0,2). The upper graph shows results for 
Wo = 0.5 and W2 = 0.5, while in the lower graph Wo = 0.01 
and W2 = 0.99 and a reweighting factor W = 99 as used 
(resulting in W^ ^ Wi ^ 0.5). 



As a first, simple example of how reweighting affects 
the statistics of a simulation we will discuss a simple 
random process. Consider a random variable n, which 
can take two different values, or 2. Let Wq and W2 
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designate the probability for these two outcomes. The 
expected fractions and standard deviations of these out- 
comes from N random selections are given by 



(^„,o) ^Wo± 



{Sua) ^W2± 



N 



N 



(34) 
(35) 



We consider an expectation value of a form similar to the 
QMC susceptibility, Eq. (p6|); 



(<5n,o) 

with a relative standard deviation 



J_ 1 W2 1 
Wo WoVWoVtv' 




(36) 



(37) 



This formula becomes valid for large N, when the stan- 
dard deviation is small. As Wo decreases the standard 
deviation increases, but one can reweight the two out- 
comes by assigning an additional weight W to the n — 
outcome such that a transition from n = to n = 2 is 
accepted with probability 1/W, while a transition from 
n = 2 to ri = is always accepted. After such a reweight- 
ing the probabilities of obtaining n = and n ~ 2 are 
given by 



and / is given by 

(/> = 



WoW 



WoW + W2 
W2 

WqW + W2 ' 



{Sn,0 + W5„ 



{Sn,o) 



(38) 
(39) 

(40) 



When calculating the standard deviation for this case 
we have to be careful since the reweighting introduces 
correlations into the system. This is clearly visualized in 
Fig. H, where in the upper graph a series of independent 
outcomes with equal probability {Wq — W2 =0.5) are 
shown , while in the lower graph a case with Wq = 0.01 
is shown with a reweighting factor of W = 99 (leading to 
W^^W!^^ 0.5). 

Let us now calculate the standard deviations for this 
case. We are interested in minimizing the standard de- 
viation of a variable evaluated in a simulation and also 
need to calculate its statistical error. In a standard MC 
simulation one usually wants to calculate the average and 
the standard deviation of the average for some quantity 
X. This is typically achieved by dividing the run into a 
number of bins, N , and saving the average of x for each 
bin. If the bins are statistically independent the final av- 
erage and standard deviation can be calculated according 
to 



1 ^ 
N 51 



N 



(41) 



and 



N 



(42) 



When studying the behavior of the standard deviation 
itself, we also want to obtain an estimate of its accu- 
racy. This can be done by dividing the N bins into M 
sets containing N/M bins each. For each set a standard 
deviation ctx can be calculated according to 



= Vx^ 



(43) 



where the bar denotes an average of the N/M bins within 
the set. The final standard deviation and its statistical 
fluctuation are then given by 



1 

i=l 



(44) 



and 



M 



(45) 



Eq. ( [l4| ) represents the standard deviation of the dis- 
tribution of the binned values a;, and not the standard 
deviation of an average of these. It does not decrease 
as the number of measurements is increased, but it is 
dependent on the number of MC steps in each bin, A^bin, 
and will decrease as l/\/Abin- Hence it is important to 
state the number of MC steps in the bins for which the 
deviation is calculated. The statistical error of this stan- 
dard deviation, acr^ , will, on the other hand, decrease as 

1/Vn. 

In this manner we can calculate the standard devia- 
tion for /. In order to show the standard deviation as a 
function of the reweighted probability Wq, Eq. (|39|), can 
be inverted to express the necessary weight factor that 
causes the average to change from Wo to Wq, 



W 



W^{l-Wo) 
Woil 



(46) 



Simulation results for the standard deviation of / as a 
function of Wq is shown in Fig. ^. We see that the 
reweighting actually increases the standard deviation. 
This is due to the rapidly increasing auto correlation 
times. The autocorrelation function 



Csit) 



((5„,o(^ + t)(5„,o(^)) - (<5„.o(z))^ 



(47) 



is shown in Fig. |l^, and one can see that the autocorrela- 
tion times (inversely proportional to the slopes in Fig. |l^) 
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FIG. 9. The relative standard deviation deviation of / 
as a function of Wq, Calculations are performed over bins 
containing Ny^i^ = 2000 MC steps. 
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FIG. 10. Autocorrelation function, Eq.(p7|), for different 
values of Wq. Shown are resuHs for Wq = 0.02, 0.05, 0.1, 0.2, 
0.3, 0.4, and 0.5 (curves left to right). 



are proportional to Wq. Notice that the longest auto- 
correlation times are significantly shorter than the indi- 
vidual bins (consisting of 2000 MC steps) used above, a 
criteria for the analysis to be valid. 

This simple example seems to indicate that reweight- 
ing does not decrease the statistical errors. However, in 
a standard Monte Carlo simulation the measured quan- 
tities are not independent even with no reweighting, and 
formula Eq. ( p7|) contains measured quantities different 
from Eq. ( |3^ ) considered above. Therefore the reweight- 
ing will affect autocorrelation times differently than in 
the above example, jand reweighting can actually decrease 
the standard error .E3 

Using the above technique we can study how the rel- 
ative error in the susceptibility of the semi-frustrated 
model is affected by reweighting. An initial run without 
reweighting has to be done first to determine the aver- 
age sign (Snfi) — Wq- Thereafter Eq. ( p6[ ) can be used 
to determine the desired weight factors. In Fig. 11 the 
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FIG. 11. The average sign as function of system volume 
at temperature T/J — 1. Shown are results of an unre- 
stricted simulation (circles), and a simulation restricted to the 
0- and 2-meron sectors without reweighting (squares). The 
line shows a slope of -2. 



average sign in a simulation of the semifrustrated model 
with J(l) = J(\/2) = J is shown as a function of lattice 
volume = L X L at a temperature T/J — 1.0. For 
comparison we first performed a standard simulation by 
sampling all meron sectors, which leads to a severe sign 
problem with a sign that decreases exponentially in sys- 
tem volume. Next we sampled only the zero- and two- 
meron sectors without reweighting, which dramatically 
increases the average sign. The scaling changes from ex- 
ponential to quadratic in the volume, as can clearly be 
seen from the graph. 
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FIG. 12. Relative standard deviation of the z-susceptibility 
of the semifrustrated model as a function of reweighting. 
Shown are results for systems of linear size 8 (circles), 12 
(squares), 16 (diamonds), and 20 (triangles). The standard 
deviation is calculated for bins containing A^bin — 10^ MC 
steps at temperature T/J — 1.0. 

Having determined the average sign without reweight- 
ing we now use Eq. ( p6[ ) to determine the desired weight 
factors. In Fig. O the standard deviation (H) of the 
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susceptibility, calculated using bins containing 1000 MC 
steps. Results are shown for systems of linear size N = 
8,12,16 and 20 at temperature T/J = 1.0. Reweight- 
ing clearly helps to reduce the standard deviation, and 
there is a definite minimum in all these curves indicating 
an optimal reweighting. The optimally reweighted sign 
always appears to be less than 0.5, and decreases with 
decreasing sign (and increasing volume) . 

Having determined that there is an optimal reweight- 
ing we will next consider whether reweighting changes the 
scaling with system size of the relative statistical error. 
Let us first consider how the standard deviation scales 
with no reweighting. Since the sign decreases quadrat- 
ically in the volume V we can derive the scaling of the 
relative error in the sign, under the ideal (and typically 
false) assumption that individual measurements are com- 
pletely independent. Using that s ~ we arrive at 



and this indicates that in this case the statistics needed 
actually increases linearly in system volume, and not 
quadratically as stated in Ref. 

0.20 I ^ . ^ . ^ . ^ . 1 
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FIG. 13. The standard deviation of the susceptibility as a 
function of system volume V. Shown are: the z-component for 
the semifrustrated model without (circles) and with (squares) 
optimal reweighting, the a;-component for the semifrustrated 
model (diamonds) and the rotationally invariant susceptibil- 
ity of the ferromagnetic model (triangles). Statistical errors 
are smaller than symbol size. 

In order to study the scaling, the standard deviation 
for bins containing A^bin = 10"' MC measurements of the 
susceptibility is shown in Fig. |l^. Four susceptibilites are 
shown: the 2-componcnt for the semifrustrated model 
without and with optimal reweighting, the x-component 
for the semifrustrated model and the rotationally invari- 
ant susceptibility of the isotropic ferromagnetic model. 
The two latter quantities can be obtained in simulations 
without sign problems, as discussed above. 

Let us first consider the z-component of the suscep- 
tibility for the semifrustrated model without and with 



optimal reweighting. For both cases the graph suggests 
a linear increase in relative error. We have to keep in 
mind that Eq. ( ^ ) does not have to be valid, since there 
are autocorrelations in the simulation, and the results 
in Fig. |l^ do not exclude that the scaling changes when 
approaching the thermodynamic limit (due to increas- 
ing autocorrelation times), but both results do support 
an approximately linear increase. It appears that the 
reweighting in this case changes only the prefactor of 
the volume scaling, but not the exponent. Whether this 
means that reweightinffi-|has completely eliminated the 
remaining sign problcmcj is not clear (the error remains 
larger than that for the ferromagnet susceptibility), but 
it is clear that the reweighting reduces the standard de- 
viations by a significant factor. In any case, an algorithm 
that changes the functional dependence of the size scal- 
ing of the statistics from exponential to polynomial can 
be considered a solution to the problem. 

The x-susceptibility of the semifrustrated model, 
which is evaluated with an algorithm without sign prob- 
lems, shows a constant standard deviation, which may be 
related to the fact that the susceptibility itself has con- 
verged to its thermodynamic limit for these system sizes 
(see next section). For the isotropic ferromagnetic model, 
the susceptibility still shows a linearly increasing error, 
but as already noted the slope is much smaller than for 
the semifrustrated case. 

This concludes our discussion of the reweighting tech- 
nique. In future work it would be interesting to ex- 
plore how the optimal reweighting can be determined 
directly from quantities measured during one single test 
run, rather than by explicitly measuring the standard 
deviations as we have done here. 



V. RESULTS 

In this section we will present results for the semifrus- 
trated and isotropic ferromagnetic models. We will 
demonstrate that it is feasible to obtain accurate re- 
sults for large systems in the 2;-basis by using the meron- 
solution. The main motivation for this study is to analyze 
the meron-solution, and we will only briefly comment on 
the physics of the semifrustrated model and how it differs 
from the isotropic ferromagnet. We will primarily con- 
sider the semifrustrated model with J(l) = J(-\/2) = J. 
To verify the correctness of our codes we have compared 
simulation results with exact diagonalization data for sys- 
tems with 4x4 spins. 

In Fig. |l| the energy is shown for the semifrustrated 
and ferromagnetic model. The statistical errorbars are 
smaller than the symbol size and the results have con- 
verged to the thermodynamic limit. The largest system 
size used had 128 x 128 spins. It can be seen that thermal 
fluctuations and finite-temperature quantum fluctuations 
more effectively destroy the ferromagnetic correlations 
for the isotropic ferromagnet than for the semifrustrated 
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FIG. 14. Energy for the semifrustrated (filled circles) and 
ferromagnetic (empty circles) model. Statistical errorbars are 
not visible in the plot, and the curves have converged to the 
thermodynamic limit. 
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FIG. 16. The x-component of the susceptibility for the 
semifrustrated model, at ratios J(\/2)/J(l) ~ 1 (empty cir- 
cles) and J{\^)/ ,7(1) — (filled circles). Statistical errorbars 
are not visible. 
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FIG. 15. The z-component of the susceptibility for 
semifrustrated (filled symbols) and ferromagnetic model 
(empty symbols). The size effects are shown for linear system 
sizes = 16 (circles), 32 (squares) and 64 (diamonds, shown 
only for the ferromagnetic model). Statistical errorbars are 
not visible on this scale. 

The z-component of susceptibility is shown for both 
models in Fig. The low-temperature susceptibility for 
finite-size systems will approach (3N/4: for the semifrus- 
trated model, and j3N/\2 for the ferromagnetic model 
(due to rotational averaging in the latter case). This 
can be clearly seen from Fig. |l^, where the susceptibility 
is multiplied by temperature, so that a straight line at 
low temperatures indicates a Curie divergence. In the 
thermodynamic limit the uniform susceptibility should 
diverge exponentially for these models, but this can not 
be seen in Fig. |l^ due to the strong finite-size effects. 
In Fig. 16 the x-coniponent of the susceptibility for 

This quantity has 



the semifrustrated model is shown, 
converged to its thermodynamic limit, and it exhibits a 



Curie divergence at low temperatures. The ground state 
value is dependent on the next-nearest neighbor coupling 
J(\/2), as is clearly demonstrated by plotting two dif- 
ferent ratios J\/2/J(l) = 1, and J\/^J(1) = (only 
nearest neighbor interactions) in Fig. 09. 



VI. CONCLUSION 

We have studied a recently introduced meron- 
solutiodHj to the sign problem within the SSE method. 
We investigated the sign problem arising in frustrated 
spin systems and showed that the meron-solution can be 
a applied to a particular model. The problems arising 
when applying the meron solution to general models of 
frustrated spins were discussed. We found that loop algo- 
rithm typically are not ergodic and merons do not exist. 
The sign problem then persists. 

For models where the meron solution works we showed 
that the sign problem can be completely eliminated for 
certain variables and largely eliminated for other vari- 
ables. This total and partial elimination comes form 
a mapping of positive- to negative-weight contributions 
and involves no approximation. For the variables where 
the sign problem is almost eliminated the statisticral er- 
rors can be reduced using a reweighting technique.E2l We 
study how the relative statistical error behaves as a func- 
tion of lattice size with and without reweighting, and 
showed that the reweighting does not change the scaling 
behavior but typically significantly reduces the over-all 
magnitude of the fluctuations. 

It is evident that the meron solution suffers from the 
same problem in both frustrated spin and fcrmionic sys- 
tems — it is confined to a few special cases. It would be 
of great importance to be able to extend it to more gen- 
eral cases, e.g., by working in basis where the required 
loop structure appears. The possibility of finding such 
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bases should be explored. 
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